#ifndef LIB_CELL_HARMONIC_H_
#define LIB_CELL_HARMONIC_H_
#include<math.h>
#include<Mesh/_Func_.h>
#include<Eigen/IterativeLinearSolvers>
#include<Eigen/Sparse>
#include<Eigen/Dense>
#include<Mesh/libcell_iterator.h>
#include<Math/LB_Math.h>

#include<iostream>

#define quote lib_cell_quote

//#include<MatOp/SparseGenMatProd.h>
//#include<Tensors/Antisymmetric_Tensor.h>
//#include<Tensors/Tensors_Operation.h>
//#include<tools/double_tools.h>

#ifndef M_PI
#define M_PI 3.1415926
#endif
typedef struct Ham_v_prop
{
    int index;
    double ux,uy;
}Ham_v_prop;
typedef struct Ham_f_prop
{
    double weight;
}Ham_f_prop;


static inline void ham_v_prop_init(Ham_v_prop* p)
{   
    p->index=-1;
    p->ux=0;p->uy=0;
}
static inline void ham_f_prop_init(Ham_f_prop* p)
{   
    p->weight=0;
}
/*static void Ham_v_prop_init(template_v*v)
{
    Ham_v_prop* hvp=(Ham_v_prop*)malloc(sizeof(Ham_v_prop));
    hvp->index=0;
    v->user_prop=hvp;
}*/
/*static void  Ham_v_prop_free(template_v*v)
{
    Ham_v_prop* hvp=(Ham_v_prop*)(v->user_prop);
    free(hvp);
}*/
static double cotan(double *p0,double* p1,double* p2)
{
    double q1[3],q2[3];
    for(int i=0;i<3;i++)
    {
            q1[i]=p1[i]-p0[i];
            q2[i]=p2[i]-p0[i];
    }

    double temp=0,re=0;
    temp=q1[1]*q2[2]-q1[2]*q2[1];
    re+=temp*temp;
    temp=-q1[0]*q2[2]+q1[2]*q2[0];
    re+=temp*temp;
    temp=q1[0]*q2[1]-q1[1]*q2[0];
    re+=temp*temp;
    re=sqrt(re);
    temp=0;
    for(int i=0;i<3;i++)
    {
        temp+=q1[i]*q2[i];
    }   

    return temp/re; 
}
double get_cotan_from_halfface(Mesh* mesh,template_hf* hf)
{
    double re=0;
    template_c* c=hf->cell;
    if(c==NULL)
    {
        return re;
    }
    template_v* vs[3];
    for(int i=0;i<c->vertices_size;i++)
    {
        bool flag=false;
        for(int j=0;j<hf->vertices_size;j++)
        {
            if(c->vertices[i]==hf->vertices[j])
            {
                    flag=true;
            }
        }
        if(!flag)
        {
            vs[0]=c->vertices[i];
            break;
        }
    }
    vs[1]=hf->vertices[0];vs[2]=hf->vertices[1];

    re=cotan(vs[0]->point,vs[1]->point,vs[2]->point); 
    return re;
}





class Harmonic
{
public:
    Mesh mesh; 
    //Eigen::SparseMatrix<double> H;
    Eigen::MatrixXd V,LLE,F;
    // double scale;
    // double center[3];

    Harmonic()
    {
        Mesh_init(&mesh);
        mesh.dimension=2;
        mesh.manifold_require=1;
        mesh.simplex=1; 

    }
    ~Harmonic()
    {
    
    }

    
    void init()
    {
        int i = 0, j = 0;
        Ham_f_prop*fp=NULL;
        Ham_v_prop*vp=NULL;
        for(auto vit=mesh.v_begin(&mesh);vit!=mesh.v_end(&mesh);vit++)
        {
            SAFE_FREE(quote(vit)->user_prop);
            
            vp=(Ham_v_prop*)malloc(sizeof(Ham_v_prop));
            ham_v_prop_init(vp);
            quote(vit)->user_prop=vp;
            if(mesh.vertex_is_boundary(&mesh,*vit))
            {
                vp->index=j;
                j++;
            }
            else {
                vp->index=i;
                i++;
            }
        }
        for(auto fit=mesh.f_begin(&mesh);fit!=mesh.f_end(&mesh);fit++)
        {
           
            SAFE_FREE(quote(fit)->user_prop);
            
            fp=(Ham_f_prop*)malloc(sizeof(Ham_f_prop));
            ham_f_prop_init(fp);
            quote(fit)->user_prop=fp;
        }
        printf("i:%d j:%d\n",i,j);
        //printf_s("%d  %d",mesh.n_vertices(),j);
        LLE.resize(i,j);
        LLE.setZero();
    }

    void init_weight()
    {
        Ham_f_prop * fp=NULL;
        for(auto fit=mesh.f_begin(&mesh);fit!=mesh.f_end(&mesh);fit++)
        {
            fp=(Ham_f_prop*)(quote(fit)->user_prop);
            fp->weight=0;
            for(int i=0;i<2;i++)
            {
                if(quote(fit)->halffaces[i]->cell!=NULL)
                {
                    fp->weight+=get_cotan_from_halfface(&mesh,quote(fit)->halffaces[i]); 
                }   
            }
        }
    }
   
    void dense2sparse(Eigen::MatrixXd &d,Eigen::SparseMatrix<double> & s)
    {
        printf("begin dense2sparse\n");
        int num=0;
        std::vector<Eigen::Triplet<double>> tripletList;
        for(int i=0;i<d.rows();i++)
        {
            for(int j=0;j<d.cols();j++)
            {
                //printf("%lf ",d.coeff(i,j));
                //if(fabs(d.coeff(i,j))>1e-8 )
                if(d.coeff(i,j)!=d.coeff(j,i)&&num<600)
                {
                    //printf("%lf %lf\n",d.coeff(i,j),d.coeff(j,i));
                    num++;
                }

                if(d.coeff(i,j)!=0)
                {

                    tripletList.push_back(Eigen::Triplet<double>(i,j,d.coeff(i,j)));
                }
            }
        }
        s.setFromTriplets(tripletList.begin(),tripletList.end());
        printf("end dense2spase %d \n",tripletList.size());
    }

    void construct_LLE()
    {
        Eigen::SparseMatrix<double,Eigen::RowMajor> H(LLE.rows(),LLE.rows());
        H.reserve(Eigen::VectorXi::Constant(LLE.rows(),30));
        //H.setZero();
        
        template_v* vs[2]={NULL,NULL};
        Eigen::MatrixXd v(H.rows(),mesh.num_v(&mesh)-H.rows());
        v.setZero();
        LLE=v;
        Ham_v_prop* vp=NULL,*vp0=NULL;Ham_f_prop*fp1=NULL;

        for(auto fit=mesh.f_begin(&mesh);fit!=mesh.f_end(&mesh);fit++)
        {
            if(!mesh.face_is_boundary(&mesh,*fit ) )
            {
                vs[0]=quote(fit)->vertices[0];
                vs[1]=quote(fit)->vertices[1];
                fp1=(Ham_f_prop*)(quote(fit)->user_prop);
                for(int i=0;i<2;i++)
                {
                    if(!mesh.vertex_is_boundary(&mesh,*(vs[i]) ))
                    {
                        vp=(Ham_v_prop*)(vs[i]->user_prop);
                        vp0=(Ham_v_prop*)(vs[(i+1)%2]->user_prop);
                        H.coeffRef(vp->index,vp->index)+=fp1->weight;
                        
                        if(!mesh.vertex_is_boundary(&mesh,* (vs[(i+1)%2]) ))
                        {
                            H.coeffRef(vp->index,vp0->index)-=fp1->weight;

                            //H.coeffRef(vp->index,vp0->index)=-fp1->weight;
                        }
                        else
                        {
                            v.coeffRef(vp->index,vp0->index)=fp1->weight;
                        }
                    }
                }   
            }
        }   

        H.makeCompressed();

        Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>llt;
        llt.compute(H);
        for(int j=0;j<v.cols();j++)
        {
            LLE.col(j)=llt.solve(v.col(j));

            //std::cout<<LLE.col(j)<<std::endl;
        } 
    }

    void natural_confromal_boudary()
    {
        int n = mesh.num_v(&mesh),int_v=LLE.rows();
        Eigen::SparseMatrix<double,Eigen::RowMajor > temp_H(2*LLE.cols(),2*n);
        temp_H.reserve(Eigen::VectorXi::Constant(2*LLE.cols(),50));
 
        Eigen::VectorXd vx(2 * LLE.cols()), x;
        // double *p1[2]={0};
        // int ini=0;
        // for(auto vit =mesh.v_begin(&mesh);vit!=mesh.v_end(&mesh);vit++)
        // {
        //     if(ini>=2)
        //     {
        //         break;
        //     }
        //     if(mesh.vertex_is_boundary(&mesh,*vit))
        //     {
        //         p1[ini]=quote(vit)->point;
        //         ini++;
        //     }
        // }
        vx.setZero();
        vx.coeffRef(0) = 0; vx.coeffRef(0 + LLE.cols()) = 0;
        vx.coeffRef(1) = 0; vx.coeffRef(1 + LLE.cols()) = 0.2;
       

        //printf("dis:%lf\n",distance_of_two_points(p1[0],p1[2],3 ) );

        //temp_H.setZero();

        // printf_s("%d  %d\r\n",n,LL.rows()+LL.cols());
        for (auto vit=mesh.v_begin(&mesh);vit!=mesh.v_end(&mesh);vit++)
        {
            Ham_v_prop* vp=(Ham_v_prop*)(quote(vit)->user_prop);
            if (mesh.vertex_is_boundary(&mesh,*vit))
            {
                if (vp->index == 0 || vp->index == 1)
                {
                   temp_H.coeffRef(vp->index, int_v + vp->index) = 1;
                   temp_H.coeffRef(LLE.cols() + vp->index, n + int_v + vp->index) = 1;  
                }
                else
                {
                    for (auto vfit=mesh.vf_begin(&mesh,*vit);vfit!=mesh.vf_end(&mesh,*vit);vfit++)
                    {
                        Ham_f_prop* fp=(Ham_f_prop*)(quote(vfit)->user_prop );

                        temp_H.coeffRef(vp->index, int_v+vp->index) += 1.0*fp->weight;
                        temp_H.coeffRef(LLE.cols()+vp->index, n+int_v + vp->index)+= 1.0*fp->weight;
                        template_v* v= (quote(vfit)->vertices[0] == quote(vit) )? (quote(vfit)->vertices[1]): (quote(vfit)->vertices[0]);
                        Ham_v_prop* vp1=(Ham_v_prop*)(v->user_prop);
                        if (mesh.vertex_is_boundary(&mesh,*v))
                        {
                            temp_H.coeffRef(vp->index, int_v + vp1->index) = -1.0*fp->weight;
                            temp_H.coeffRef(LLE.cols()+vp->index, n+int_v + vp1->index) = -1.0*fp->weight;
                            //*****************
                            if( mesh.face_is_boundary(&mesh,*vfit) )
                            {
                                //temp_H.coeffRef(vp->index,n+int_v+vp1->index)=1;
                                if(quote(vfit)->vertices[0]==quote(vit))
                                {
                                    temp_H.coeffRef(vp->index,n+int_v+vp1->index)=1;
                                    temp_H.coeffRef(LLE.cols()+vp->index,int_v+vp1->index)=-1;

                                    //printf("1 ");
                                }
                                else if(quote(vfit)->vertices[1]==quote(vit))
                                {
                                    temp_H.coeffRef(vp->index,n+int_v+vp1->index)=-1;
                                    temp_H.coeffRef(LLE.cols()+vp->index,int_v+vp1->index)=1;

                                    //printf("2 ");

                                }

                            }
                            //**************************
                        }
                        else
                        {
                            temp_H.coeffRef(vp->index,  vp1->index) = -1.0*fp->weight;
                            temp_H.coeffRef(LLE.cols() + vp->index, n + vp1->index) = -1.0*fp->weight;
                        }                    
                    }
                    //printf("\n");
                }
            }
        }

        Eigen::MatrixXd temp_M(2*LLE.cols(),2*LLE.cols());
        temp_M.setZero();

        //Eigen::SparseMatrix<double > temp_M(2*LLE.cols(),2*LLE.cols());
       // printf("%d %d %d %d\n",temp_H.rows(),temp_H.cols(),temp_L.rows(),temp_L.cols());
        //Eigen:: SparseMatrix<double > temp_M1=temp_H.block(0,0,LLE.cols(),LLE.rows())*LLE;

        temp_M.block(0,0,LLE.cols(),LLE.cols())=temp_H.block(0,0,LLE.cols(),LLE.rows())*LLE+temp_H.block(0,LLE.rows(),LLE.cols(),LLE.cols());

        temp_M.block(0,LLE.cols(),LLE.cols(),LLE.cols())=temp_H.block(0,n,LLE.cols(),LLE.rows())*LLE+temp_H.block(0,n+LLE.rows(),LLE.cols(),LLE.cols());

        temp_M.block(LLE.cols(),0,LLE.cols(),LLE.cols())=temp_H.block(LLE.cols(),0,LLE.cols(),LLE.rows())*LLE+
            temp_H.block(LLE.cols(),LLE.rows(),LLE.cols(),LLE.cols());

        temp_M.block(LLE.cols(),LLE.cols(),LLE.cols(),LLE.cols())=temp_H.block(LLE.cols(),n,LLE.cols(),LLE.rows())*LLE+
            temp_H.block(LLE.cols(),n+LLE.rows(),LLE.cols(),LLE.cols());

        printf("begin natural solve %d %d\n",LLE.cols(),LLE.rows());
        ///std::cout<<LLE<<std::endl;

        
        //Eigen::SparseMatrix<double >temp_MM(2 * LLE.cols(), 2 * LLE.cols());

        //dense2sparse(temp_M,temp_MM);

        // Eigen::SparseLU<Eigen::SparseMatrix<double>> LU;
        // // //Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>llt;
        // LU.compute(temp_MM);
        // x = LU.solve(vx);
        
        Eigen::PartialPivLU<Eigen::MatrixXd> luSolver;
        luSolver.compute(temp_M);
        x=luSolver.solve(vx);
        
        //std::cout<<x<<std::endl;
        //Eigen::BiCGSTAB<Eigen::SparseMatrix<double>, Eigen::IncompleteLUT<double>> solver;
        // solver.compute(temp_MM);
        // x = solver.solve(vx);
        // //x = temp_M.lu().solve(vx);
        // printf_s("%lf\r\n", vx.coeffRef(1 + LL.cols()));
        for(auto vit= mesh.v_begin(&mesh);vit!=mesh.v_end(&mesh);vit++)
        {
            if(mesh.vertex_is_boundary(&mesh,*vit))
            {
                Ham_v_prop* vp=(Ham_v_prop*)(quote(vit)->user_prop);
                vp->ux=x.coeff(vp->index);
                vp->uy=x.coeff(LLE.cols()+vp->index);
            }
        }
      
        printf("end mat hrerere\n");

    }

    void map_squre()
    {
        mesh.external_cell_init_(&mesh);

        int size=node_size(mesh.external_cell.halffaces);
        int step=size/4;
        double dir[4][2];
        dir[0][0]=1;dir[0][1]=0;    
        dir[1][0]=0;dir[1][1]=1;    
        dir[2][0]=-1;dir[2][1]=0;   
        dir[3][0]=0;dir[3][1]=-1;
        int i=1,k=0;

    //  printf("%d\n",step);
        template_v* v0=NULL,*v1=NULL;
        Ham_v_prop* temp=NULL,*temp1=NULL;
        int steps[4];
        steps[0]=step;steps[1]=2*step;steps[2]=3*step;steps[3]=size;
        for(Node* nit=mesh.external_cell.halffaces;nit!=NULL;nit=(Node*)(nit->Next))
        {
            v0=((template_hf*)(nit->value))->vertices[1];
            v1=((template_hf*)(nit->value))->vertices[0];
            //v=(template_v*)(nit->value);
            temp=(Ham_v_prop*)(v0->user_prop);
            temp1=(Ham_v_prop*)(v1->user_prop);
            
            if(k==3)
            {
                double temp_step=size-step*3;
                temp1->ux=temp->ux+dir[k][0]/(double)temp_step;
                temp1->uy=temp->uy+dir[k][1]/(double)temp_step;
            }
            else
            {
                temp1->ux=temp->ux+dir[k][0]/(double)step;

                temp1->uy=temp->uy+dir[k][1]/(double)step;
            }
        
            //printf("v1:%d v2:%d,ux:%lf uy:%lf\n",v0->id,v1->id,temp1->ux,temp1->uy);
            if(i==steps[0]||i==steps[1]||i==steps[2]||i==steps[3])
            {
                k++;
            }
            i++;
        }   
        temp1->ux=0;temp1->uy=0;
    }

    void map_disk()
    {
        int i=0;
        double *dis=(double*)malloc(sizeof(double)*LLE.cols());
        memset(dis,0,sizeof(double)*LLE.cols());

        mesh.external_cell_init_(&mesh);
        for(auto chfit= mesh.chf_begin(&mesh,mesh.external_cell);
            chfit!=mesh.chf_end(&mesh,mesh.external_cell);chfit++ )
        {
            dis[i]= distance_of_two_points(quote(chfit)->vertices[0]->point,
                quote(chfit)->vertices[1]->point,3)+ dis[(i-1+LLE.cols())%LLE.cols() ];
            i++;
        }
        double sum = dis[i-1];
        i = 0; double theta = 0;

        for(auto chfit= mesh.chf_begin(&mesh,mesh.external_cell);
            chfit!=mesh.chf_end(&mesh,mesh.external_cell);chfit++ )
        {
            theta = 2*M_PI*dis[i]/sum;
            Ham_v_prop* vp=(Ham_v_prop*) (quote(chfit)->vertices[0]->user_prop);
            vp->ux=cos(theta);
            vp->uy=sin(theta);
            i++;
        }
    }
    
    void map()
    {
        Eigen::VectorXd ux(LLE.cols()),uy(LLE.cols()),x(LLE.cols()),y(LLE.cols());
        Ham_v_prop* temp=NULL;
        for(auto iter=mesh.v_begin(&mesh);iter!=mesh.v_end(&mesh);iter++)
        {
            if(mesh.vertex_is_boundary(&mesh,*(iter)))
            {
                temp=(Ham_v_prop*)(quote(iter)->user_prop);
                ux.coeffRef(temp->index)=temp->ux;
                uy.coeffRef(temp->index)=temp->uy;

            }
        }   
        x=LLE*ux;y=LLE*uy;

        for(auto iter=mesh.v_begin(&mesh);iter!=mesh.v_end(&mesh);iter++)
        {
            if(!mesh.vertex_is_boundary(&mesh,*(iter)))
            {
                temp=(Ham_v_prop*)(quote(iter)->user_prop);
                temp->ux=x.coeff(temp->index);
                //printf("ux:%lf\n",temp->ux);
                temp->uy=y.coeff(temp->index);

            }
        }   
    }

   
};

class Geometry_image
{
public:
    Geometry_image(Mesh *m):mesh(m)
    {
        image_tensor=NULL;
        test_data=NULL;
        image_tensor1=NULL;
        image_tensor2=NULL;
    }
    ~Geometry_image()
    {
        if(image_tensor!=NULL)
        {
            free(image_tensor);
        }
        if(test_data!=NULL)
        {
            free(test_data);
        }
        if(image_tensor1!=NULL)
        {
            free(image_tensor1);
        }
        if(image_tensor2!=NULL)
        {
            free(image_tensor2);
        }
    }
    void WriteBMP(unsigned char*img,const char* filename,int w,int h)
    {
#ifdef WIN32
    FILE *fp;
    fopen_s(&fp,filename,"r");
#else
    FILE *fp=fopen(filename,"wb");
#endif
        int l=(w*3+3)/4*4;
        int bmi[]= {l*h+54,0,54,40,w,h,1|3*8<<16,0,l*h,0,0,100,0};
        //FILE *fp = fopen(filename,"wb");
        fprintf(fp,"BM");
        fwrite(&bmi,52,1,fp);
        fwrite(img,1,l*h,fp);
        fclose(fp);
    }
    void double_image(int w,int h)
    {
        image_tensor1=(unsigned char*)malloc(sizeof(unsigned char)*w*h*3);
        image_tensor2=(unsigned char*)malloc(sizeof(unsigned char)*w*h*3);  
        for(int i=0;i<w;i++)
        {
            for(int j=0;j<h;j++)
            {
                for(int k=0;k<3;k++)
                {
                    int temp=round(test_data[(j*w+i)*3+k]*(255*255-1));
                    int w1=temp/255;
                    int w2=temp-w1*255;
                    image_tensor1[(j*w+i)*3+k]=w1;
                    image_tensor2[(j*w+i)*3+k]=w2;

                }
            }
        }
    }

    void double_image2mesh(unsigned char *data1,unsigned char *data2,int w,int h,Mesh* m)
    {
        Mesh_init(m);
        m->simplex=1;
#ifdef MANIFOLD_REQUIRE
        m->dimension=2;
#endif
        double p[3];
        printf("hereli %d %d \n",w,h);
        template_v**indexs=(template_v**)malloc(sizeof(template_v*)*w*h);
        printf("hereli\n");
        for(int i=0;i<w;i++)
        {
            for(int j=0;j<h;j++)
            {
                for(int k=0;k<3;k++)
                {
                    p[k]=(data1[(j*w+i)*3+k]*255+data2[(j*w+i)*3+k])/(double)(255.0*255.0-1.0);
                    //p[k]=test_data[(j*w+i)*3+k];
                }
                template_v*v=m->create_vertexv(m,p,3);
                indexs[j*w+i]=v;    
            }
        }
        template_v* vs[3];
        for(int i=0;i<w-1;i++)
        {
            for(int j=0;j<h-1;j++)
            {
                vs[0]=indexs[j*w+i];
                vs[2]=indexs[(j+1)*w+i];
                vs[1]=indexs[(j+1)*w+i+1];

                m->create_cellv(m,vs,3);
                vs[0]=indexs[j*w+i];
                vs[2]=indexs[(j+1)*w+i+1];
                vs[1]=indexs[j*w+i+1];
                m->create_cellv(m,vs,3);    
            }
        }
        free(indexs);   
    }
    void  image2mesh(unsigned char * data,int w,int h,Mesh* m)
    {
        Mesh_init(m);
        m->simplex=1;
#ifdef MANIFOLD_REQUIRE
        m->dimension=2;
#endif
        double p[3];
        printf("hereli %d %d \n",w,h);
        template_v**indexs=(template_v**)malloc(sizeof(template_v*)*w*h);
        printf("hereli\n");
        for(int i=0;i<w;i++)
        {
            for(int j=0;j<h;j++)
            {
                for(int k=0;k<3;k++)
                {
                 
                    p[k]=data[(j*w+i)*3+k]/255.0;
                }
                template_v*v=m->create_vertexv(m,p,3);
                indexs[j*w+i]=v;    
            }
        }
        template_v* vs[3];
        for(int i=0;i<w-1;i++)
        {
            for(int j=0;j<h-1;j++)
            {
                vs[0]=indexs[j*w+i];
                vs[2]=indexs[(j+1)*w+i];
                vs[1]=indexs[(j+1)*w+i+1];

                m->create_cellv(m,vs,3);
                vs[0]=indexs[j*w+i];
                vs[2]=indexs[(j+1)*w+i+1];
                vs[1]=indexs[j*w+i+1];
                m->create_cellv(m,vs,3);    
            }
        }
        free(indexs);   

    }
    void regularize_uv()
    {
        double max_x=0,min_x=1,max_y=0,min_y=1;
        Ham_v_prop* temp=NULL;      
        for(auto it=mesh->v_begin(mesh);it!=mesh->v_end(mesh);it++)
        {

                temp=(Ham_v_prop*)(quote(it)->user_prop);
                if(temp->uy>max_y)
                {
                    max_y=temp->uy;
                }
                if(temp->uy<min_y)
                {
                    min_y=temp->uy;
                }
                if(temp->ux>max_x)
                {
                    max_x=temp->ux;
                }
                if(temp->ux<min_x)
                {
                    min_x=temp->ux;
                }
        }
        printf("%lf %lf %lf %lf\n",max_x,min_x,max_y,min_y);    

    }
    //边界调整
    void map_squre(int w,int h)
    {
        mesh->external_cell_init_(mesh);
        template_v* v0=NULL,*v1=NULL,*vv1=NULL,*vv2=NULL;   
        Ham_v_prop* temp=NULL,*temp1=NULL;
        int x1,x2,y1,y2;    
        for(Node* nit=mesh->external_cell.halffaces;nit!=NULL;nit=(Node*)(nit->Next))
        {
            v0=((template_hf*)(nit->value))->vertices[1];
            v1=((template_hf*)(nit->value))->vertices[0];
            //v=(template_v*)(nit->value);
            temp=(Ham_v_prop*)(v0->user_prop);
            temp1=(Ham_v_prop*)(v1->user_prop);
            x1=round(temp->ux*(w-1));x2=round(temp1->ux*(w-1));
            y1=round(temp->uy*(h-1));y2=round(temp1->uy*(h-1));
            double loc=0;
            double lambda=0;
            if(x1!=x2)
            {
                vv1=v0;vv2=v1;
                int in1=x1,in2=x2;
                if(x2>x1)
                {

                }
                else
                {
                    vv1=v1;vv2=v0;
                    in1=x2,in2=x1;
                }

                for(int i=in1;i<=in2&&i<w;i++)
                {
                    temp=(Ham_v_prop*)(vv1->user_prop);temp1=(Ham_v_prop*)(vv2->user_prop);

                    lambda=(i/(double)(w-1)-temp->ux)/(temp1->ux-temp->ux);
                    if(lambda>=0&&lambda<=1)
                    {

                        for(int k=0;k<3;k++)
                        {
                            loc=(lambda*vv2->point[k]+(1-lambda)*vv1->point[k]);
                            int temp_i=round(loc*254.0);
                    
                            test_data[(y1*w+i)*3+k]=loc;
                            image_tensor[(y1*w+i)*3+k]=temp_i;
                        }
                    }   
                    //image_tensor[(y1*w+i)*3];
                    //printf("\n");
                }
            }
            if(y1!=y2)
            {
                vv1=v0;vv2=v1;
                int in1=y1,in2=y2;
                if(y2>y1)
                {

                }
                else
                {
                    vv1=v1;vv2=v0;
                    in1=y2,in2=y1;
                }
                for(int i=in1;i<=in2&&i<h;i++)
                {
                    temp=(Ham_v_prop*)(vv1->user_prop);temp1=(Ham_v_prop*)(vv2->user_prop);

                    lambda=(i/(double)(h-1)-temp->uy)/(temp1->uy-temp->uy);
                    if(lambda>=0&&lambda<=1)
                    {


                        for(int k=0;k<3;k++)
                        {
                            loc=(lambda*vv2->point[k]+(1-lambda)*vv1->point[k]);
                            int temp_i=round(loc*254.0);
                 
                            test_data[(i*w+x1)*3+k]=loc;
                            image_tensor[(i*w+x1)*3+k]=temp_i;
                        }   
                    }
                    //printf("\n");
                }

            }
            
        }   

    }
    void regular_sample(int w,int h)
    {
        regularize_uv();
        image_tensor=(unsigned char*)malloc(sizeof(unsigned char)*w*h*3);
        test_data=(float*)malloc(sizeof(float)*w*h*3);
        memset(test_data,0,sizeof(float)*w*h*3);
        memset(image_tensor,0,sizeof(unsigned char)*w*h*3);
        double max_y=0,min_y=1;
        double max_x=0,min_x=1;
        Eigen::MatrixXd A(3,3);
        A.coeffRef(2,0)=1;A.coeffRef(2,1)=1;A.coeffRef(2,2)=1;  
        Eigen::VectorXd b(3);
        b.coeffRef(2)=1;
        Ham_v_prop* temp=NULL;
        int max_y_,min_y_,max_x_,min_x_;
        for(auto cit=mesh->c_begin(mesh);cit!=mesh->c_end(mesh);cit++)
        {
            //printf("begin\n");
            max_y=0;min_y=1;
            max_x=0;min_x=1;
            int max_y_index=0,min_y_index=0;
            int max_x_index=0,min_x_index=0;
            for(int j=0;j<quote(cit)->vertices_size;j++)
            {
                temp=(Ham_v_prop*)(quote(cit)->vertices[j]->user_prop);
                if(temp->uy>max_y)
                {
                    max_y_index=j;
                    max_y=temp->uy;
                }
                if(temp->uy<min_y)
                {
                    min_y_index=j;
                    min_y=temp->uy;
                }
                if(temp->ux>max_x)
                {
                    max_x_index=j;
                    max_x=temp->ux;
                }
                if(temp->ux<min_x)
                {
                    min_x_index=j;
                    min_x=temp->ux;
                }
                A.coeffRef(0,j)=temp->ux;
                A.coeffRef(1,j)=temp->uy;
            }
            min_y_=((int)(min_y*h));
            max_y_=((int)(max_y*h));    
            min_x_=((int)(min_x*w));
            max_x_=((int)(max_x*w));    

            if(min_y_==max_y_)
            {
                for(int i=0;i<3;i++)
                {
                    A.coeffRef(1,i)=0;
                    if(i!=min_x_index&&i!=max_x_index)
                    {
                        A.coeffRef(1,i)=1;
                    }
                }
                //A.coeffRef(1,0)=1;A.coeffRef(1,1)=0;A.coeffRef(1,2)=0;
            }   
            else if(min_x_==max_x_)
            {
                //b.coeff(1)=0;
                for(int i=0;i<3;i++)
                {
                    A.coeffRef(0,i)=0;
                    if(i!=min_y_index&&i!=max_y_index)
                    {
                        A.coeffRef(0,i)=1;
                    }
                }
            }
            {
                Eigen::MatrixXd Ain=A.inverse(),x,y;
                //printf("%d %d %d %d\n",min_y_,max_y_,min_x_,max_x_);      
                for(int i=min_x_;i<=max_x_&&i<w&&i>=0;i++)
                {

                    for(int j=min_y_;j<=max_y_&&j<h&&j>=0;j++)
                    {
                    //  printf("libo\n");
                        b.coeffRef(0)=i/(double)(w-1);
                        b.coeffRef(1)=j/(double)(h-1);
                        //bool flag= false;
                        if(min_y_==max_y_)
                        {
                        
                            b.coeffRef(1)=0;
                            //x=Ain*b;
                            //y=A*x;
                        //  flag=true;
                        }
                        else if(min_x_==max_x_)
                        {
                            b.coeffRef(0)=0;
                        //  flag=true;
                        }
                        else
                        {

                        }
                    //  printf("b %lf %lf\n",b.coeff(0),b.coeff(1));    
                        x=Ain*b;
                        y=A*x;
                    //  printf("i:%d j:%d jiance x:%lf %lf %lf b:%lf %lf,A:%lf %lf %lf %lf\n",i,j,x.coeff(0),x.coeff(1),x.coeff(2),b.coeff(0),b.coeff(1),A.coeff(1,0),A.coeff(1,1),A.coeff(1,2),A.coeff(1,0));
                        if(x.coeff(0)>=-0.00001&&x.coeff(1)>=-0.00001&&x.coeff(2)>=-0.00001)
                        {
                    //      printf("here\n");
                            double loc[3]={0,0,0};
                            for(int k=0;k<3;k++)
                            {
                                loc[0]+=x.coeff(k)*quote(cit)->vertices[k]->point[0];
                                loc[1]+=x.coeff(k)*quote(cit)->vertices[k]->point[1];  
                                loc[2]+=x.coeff(k)*quote(cit)->vertices[k]->point[2];
                            }
                            //printf("here\n");
                            for(int k=0;k<3;k++)
                            {
                                int temp_i=round(loc[k]*254);

            //                  printf("i:%d j%d libo%d\n",i,j,(j*w+i)*3+k);
                                if((j*w+i)*3+k>=30000)
                                {
                                    printf("hauile\n");
                                    return;
                                }
                                test_data[(j*w+i)*3+k]=loc[k];
                                //printf("%lf %d ",loc[k],temp_i);
                                image_tensor[(j*w+i)*3+k]=temp_i;
                                //printf("%d ",image_tensor[(i*w+j)*3+k]);  
                            }
                //          printf("\n");
                        }                           
                    }
                
                }
                //printf("end\n");
            }
        }
    
            
        map_squre(w,h);
        double_image(w,h);  

    }

    float *test_data;
    
//protected:
    unsigned char *image_tensor;
    unsigned char *image_tensor1,*image_tensor2;
    Mesh* mesh;
};

#endif
